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Abstract 



We present a study of sediment transport in the creeping and saltation regime. 
In our model, a bed of particles is simulated with the conventional event-driven 
method. The particles are considered as hard disks in a 2d domain, with periodic 
boundary conditions in horizontal direction. The flow of the fluid over this bed 
of particles is modeled by imposing a force on each particle that depends on the 
velocity of the fluid and its height above the bed. We considered two velocity profiles 
for the fluid, parabolic and logarithmic. The first one models laminar flow and the 
second corresponds to turbulent flow. For each case we investigated the behavior 
of the saturated flux. We found that for the logarithmic profile, the saturated flux 
shows a quadratic increase with the strength of the flow, and for parabolic profile, a 
cubic increase. The velocity distribution functions are used to interpret the results. 
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1 Introduction 



The study of transport of granular material by a fluid is important for industrial 
processes as well as understanding of natural phenomena. Modeling of the sediment 
transport in rivers, as well as modeling of sand drift in the formation of dunes, 
benefits from this study. 

Saltation, surface creep, and suspension are three modes which occur during 
transportation of granular material by a fluid [I]. When the shear velocity of the 
fluid flowing over a bed of grains exceeds the friction threshold velocity for sand 
transport, the grains are driven by the fluid. At first they begin to move while still 
in continual contact with each other, yet it could happen that every now and then 
due to collisions, some particles jump by a distance of order of their diameter. This 
regime is called surface creep or reptation. As the fluid shear velocity increases, the 
particles can follow paths that take them to a height much larger than the grain 
diameter; this regime is called saltation. The grains in saltation have been named 
saltons, and the grains in creeping motion have been named reptons [5^. Suspension 
occurs at very high shear velocities, when a considerable fraction of the particles 
are transported upwards by turbulent eddies. In this regime, the grains move in the 
fiuid for long periods of time, hardly colliding with the bed or each other. Except 
for dust storms in which suspension is dominant, creeping and saltation of sands 
usually play the key role in dune formation [5] . In many of the experimental studies 
of sediment transport, grains are transported by air [3], [4], but few experiments in 
water also exist [5], [6]. 

To gain insight into the problem of sediment transport, it is important to under- 
stand the relation between the fiux of the grains transported by the fiuid and the 
velocity profile of the fluid that moves over the grains. In most instances of sediment 
transport, the flux eventually saturates at a certain strength or amplitude, u*, of 
the velocity proflle, v{y) = u^,f{y). Here, / is a function of height, y. However, 
there are situations, as at the foot of a sand dune, where the sand flux may never 
reach saturation. In such cases, where there is no saturation the variation of flux 
can be studied |T . 

Bagnold yy was first to introduce a simple fiux law, a cubic relation, expressing 
the behavior of sand fiux with the shear velocity. Other forms for the fiux law have 
also been proposed by different authors for different physical conditions [8], [9]. 

In this work, we study two profiles for the fiuid fiow, and concentrate on the 
saltation regime of the grain motion. We emphasize that for the purposes of this 
study, saltation has a different meaning than its standard usage. For instance, in 
a wind tunnel, a grain is in saltation if its trajectory is at least about 300 grain 
diameters high and at least 1000 grain diameters across. These are much larger 
than the size of the system considered here. Yet, we use this term to distinguish 
the motion from the situation where the grains constantly touch each other as they 
slowly move. Saltation is then used to mean a motion where the grains jump and 
follow a trajectory, albeit smaller than mentioned above. 

The structure of the paper is as follows. In Section II we introduce the model. 
Then in Section HI, we study the behavior of the flux as a function of the velocity 
profile, as well as the velocity distribution functions for the grains. Section IV 
is devoted to our discussions which mainly rest upon the comparison of velocity 
distributions. Finally, we present our conclusions as Section V. 

2 Simulation Model 

We use the inelastic hard sphere model lOJ. Grains are contained in a two di- 
mensional rectangular domain with periodic boundary conditions in the horizontal 
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Figure 1: The schematic representation of the system showing the particles in their 
initial state, and the velocity field at different positions in the y-direction. 

directions. The fluid flows over the bed of grains, so grains can be entrained by the 
fluid. The fluid is modeled by its velocity profile. 

The grains, modeled as disks, move under the influence of both gravity and a 
drag force that is exerted on them by the fluid. To avoid crystallizations, 20% of 
the particles have a diameter equal to 0.61q and the rest have a diameter equal to 
0.5/0 where Iq is the unit of length used throughout this paper. A gravitational 
acceleration of 12/o/to applied to all the particles, where to is the unit of time 

m- 

All the particles have the same mass, and the effect of particles on the fluid is 
neglected. 

In all cases studied, the system starts out having six layers of particles resting 
fairly compactly on each other. The fluid stands at a height that is equivalent 
to thirty two layers, 16Zo- This height is the maximum attainable by the grains; 
because in our modellization we implement a reflecting boundary at the top so that 
the particles that touch it just reverse their vertical velocity component Fig.([T]). 

2.1 Particle Motion 

We use event-driven MD [T^] to calculate the motion of the particles. As the 
particles are hard spheres, collisions take infinitesimal time and involve only two 
particles. Conservation of momentum leads to 

1 + r f ~ \- 

Ml,2 = Vi^2 T 2 yK ■ [Vl ~ V2) j K 

T i^(^(^;l-^/2))^, (1) 

where u indicates the velocities after the collision, and v denotes the velocities before 
the collision. The geometry of the collision is described by K, a unit vector pointing 
along the line connecting the centers of particle 1 to particle 2, and t is the unit 
vector in tangential direction. The energy dissipation is measured by r, the normal 
restitution coefficient, and (3, the tangential restitution coefficient. If r = I and 
/3 = ±1, collisions conserve energy and are said to be elastic. For < r < 1 or 
— 1 < /? < 1 energy is dissipated and the collisions are inelastic. We assume that 
the particles neither rotate nor roll. This is somewhat in accord with the fact that 
the sand grains are not round which makes the rolling difficult. In our simulations, 
r = 0.4 and (3 = -I. 

In using the event-driven method there are two problems in setting up the bed of 
grains: inelastic collapse, and creating a rough surface on the bottom. All particles 
after some collisions lose their energy and accumulate on the bottom and make a 
dense network of grains. So the number of collisions per unit time will diverge at 
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finite time; that is, inelastic collapse [T3j occurs. Because of the finite precision of 
the computer, multiparticle collisions can occur. For handling the inelastic collapse 
we use the Tc model with = lO^^to [TB] . 

In order to create a rough surface, r and (3 are adjusted for collisions between the 
grains and the surface. We suppose that when a particle of the first layer bounces 
against the bottom plate both the tangential and normal components of its velocity 
are reversed, i.e, r = 1 and [3 — \in Eq.([T]). In this way, the first layer is nearly fixed 
and acts as a rough surface over which other particles can move. The roughness is 
of the order of the particle diameter. 



2.2 The Effect of Fluid on the Grains 

The drag force is proportional to the difference between the particle velocity Up and 
the fiuid velocity Uf. 

F = ^{uf - Up), (2) 

where 7 is a parameter that depends on the characteristics of both the fluid and the 
grains. In laminar flow with small Reynolds number, 7 = 3mjdp, in which rj is the 
viscosity of the fiuid and dp is the diameter of the particle. We suppose that 7 = 1 
for laminar fiow. However in turbulent fiow, the fluid drag varies as the square of 
the grain speed and 7 can be written as: 

7 = \uf-Up I, (3) 

which corresponds to the Newtonian drag force where, Cd is taken from empirical 
relations, pf and pp are the density of the fluid and particle respectively. 

In general, the drag force acts on upper layers of the bed of particles and drops 
to zero for lower layers. The details depend on the velocity profile considered. 
Here, we study the dynamics of the grains for two velocity profiles: logarithmic and 
parabolic. For the parabolic profile, which models laminar flow in open channels as 
function of height, we have: 

Uf = u,{yoiyo/2-h)-y{y/2-h)), (4) 

where h is the height of fiuid in channel. This equation is written so that it satisfies 
the two boundary conditions, duf / dy — Q dX y = h, and Uf — ai y = y^. Here, y 
is the height, and j/o is the height below which the effect of the fluid on the grains 
is negligible. 

In turbulent flow, the velocity proflle of the fluid near the boundary is observed 
to be logarithmic il4j, and described by 

M/ = y'o5(y/yo), (5) 

where k is the von Karman constant. Although this relationship has been derived 
only for the region in which the shear stress is approximately constant, experiments 
show that the agreement persists through almost all of the boundary layer. In our 
simulations, yo is about the diameter of grains, yo — O.S/q. We consider no vertical 
component to the fiuid velocity. 



3 Granular Flux and Velocity Probability Distri- 
bution Function 

First we investigate the behavior of the granular flux with respect to time. Granular 
flux is the number of grains that cross the unit surface (unit line in 2d ) per time 
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Figure 2: Flux as a function of time for different values of for a logarithmic 
profile. 

unit. So the flux integrated over a cross section of the periodic domain (over height 
in 2d) is: 



where L is the length of the box and N is the total number of particles, and Ui 
denotes the horizontal component of the velocity of the ith particle. 

After some time the flux fluctuates around a constant, steady value, and sedi- 
ment transport reaches steady state. This steady state is the saturated flux, and it 
depends on the shear velocity, denoted by u*. 

For calculating the saturated flux value, we average over the flux only after the 
steady state has been reached. We expect to have a threshold velocity ut, below 
which sediment transport cannot happen. 

3.1 Logarithmic Profile 

In the aeolian case the logarithmic velocity profile is more realistic than the parabolic 
profile. When the wind blows over a rough surface, its velocity within the boundary 
layer increases logarithmically with height, Eq.([5|). The logarithmic profile has also 
been observed experimentally, for water moving over a rough surface in a channel 



Fig. ([2]) shows the variation of flux with time for different values of shear velocity, 
M,. This figure shows that transport process reaches steady state and the particle 
flux saturates after some transient time. 

The simulations could be made into movies of the grain motion. This was a 
particularly useful way of interpreting the results. In this way we estimate the 
grain motion to start at a threshold velocity of about u, = 1. With increasing 
from Ui = 1, the top layer particles begin to roll in their own layer or jump to a 
height about their diameter. This situation continues until ~ 15. This means 
that for M, < 15 most of the particles except those in the bottom layer are in the 
creeping regime. After > 15 the shear stress is enough to make some of the 
particles in the upper layers enter the saltation regime. With increasing u», the 
number of saltating particles increases, resulting in an increase in the number of 
collisions between saltons. We estimate u* ~ 15 as marking the onset of saltation in 
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Figure 3: Mean flux vs. for a logaritlimic profile, and the fit to it, Eq. pOl) 
and Bagnold's expression. Fit the Eq. (jTU]) gives a = 0.2, b = 43.9, and ut = 1.09. 
According to the simulations we estimate ttt = 1. 



this system. When > 30 enough particles have so high an energy that they move 
with the fluid stream above the other particles and have few collisions with each 
other or the rest of the grains. In this case the length of their trajectory becomes 
comparable to the system size, hence we only considered < 30. We wish to 
emphasize again that we are using the term saltation in a restricted sense in this 
study. 

The main objective of this study is to relate the saturated flux and the shear 
velocity in the regime of saltation. Our results for the behavior of the mean flux of 
particles with u^, is shown in Fig. ([3]). 

Theoretically, perhaps the most important description for saturated flux is due 
to Bagnold [T]. He found that the saturated flux at large shear velocities is given 
by: 

q^^ul (7) 
9 

Bagnold considered a mean trajectory for each grain, and supposed that the 
ejection velocity of grains from the bed scales with the shear velocity u^. This 
hypothesis is valid if the shear velocity is large enough. For small shear velocities, 
Ungar & Haff 17J supposed that height and length of the trajectory of grains is of 
the same order as the grain size, and predicted that: 

q(x pairiul-uf)^, (8) 

where d is the grain diameter. Almeida et a/.|15| also found numerically a quadratic 
description for the flux near the threshold shear velocity, 

q (X {u^ - Utf (9) 

They simulated the saltation inside a two-dimensional channel with a mobile top 
wall. Their model solves the turbulent wind fleld including feedback of the dragged 
particles. 

Our results show a slightly different quadratic dependence: 

q — a{u<, — Ut)(ua. + b) (10) 
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Figure 4: The behavior of grain velocity probabihty function for different values of 
M* for the logarithmic profile. 

where a — 0.2, 6 = 43.9, and ut — 1.09 are obtained from fitting the Eq. (fTO|) to 
our data. This value for the threshold velocity is in good agreement with what was 
estimated from the simulations; that is ut = 1. Fig. ([31) shows the fit, Eq. pI7)) . to 
the data. This indicates that both in the creeping regime and the saltating regime, 
the same quadratic function describes the behavior of flux reasonably well. Eq. (fTO|l 
predicts a stronger dependence on shear velocity, compared to Eq. ([9|) . One reason 
may be our neglect of the back action of the grains on the fluid. This effect is more 
prominent at small heights above the bed where the particle velocity differs much 
from the fluid velocity [JLSj . 

Another way to estimate the onset of saltation is to investigate the behavior 
of the grain velocity probability distribution function. The distribution function 
allows decomposition of the flux into a part due to saltation and another part due 
to creep. This interpretation is based on the fact that the reptons are slower than 
the saltons. Fig.Q shows the behavior of the grain velocity distribution for different 
values of u*. 

For small values of there is a large peak at small velocities that shows that 
all of the particles are in creeping motion. With increasing u^, the velocity of the 
moving particles increases and some of the particles enter the saltation regime. The 
uniform distribution function at higher velocities represents the saltating particles. 
At = 15 the particles can go to the saltation regime; the velocity distribution 
develops a minimum at i; ~ 1. At 20 all of particles have saltating motion. 

3.2 Parabolic Profile 

At low Reynolds numbers when a fluid flows in an open channel , its velocity varies 
with height parabolically. Now, we consider a velocity profile as in Eq.(l4|), that is 
zero at the bottom and increases parabolically with height. 

Fig. ([5]) shows the behavior of flux as a function of time. Similar to the logarith- 
mic profile, the flux reaches steady state, only now the transient time is longer. In 
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Figure 5: Flux as a function of time a for parabolic profile. 

the parabolic profile, the shear stress above the bed is greater than the logarithmic 
profile and increases rapidly with u,, so compared to the logarithmic case, the range 
of that defines the creep motion is much narrower. 

The motion starts out at about ut — 0.08 [f6 ; this is an estimate obtained from 
our simulations, and in Eq. (|12[) below is used as an input data. When = 1 the 
system is out of the creeping regime, and starts saltating. The surface y = h is 
as before a reflecting surface, so any particle that reaches it bounces back into the 
fluid. This starts at about = 1.2. 

Beyond = 2.5 this boundary produces a kind of population inversion and the 
collision of particles with each other effectively reduces the rate of increase of flux 
Fig. il). 

In fig.© we show the behavior of the mean flux with (near ut). For u» < 1 
we fitted the data with two equations below: 

q = a'{u^,-utf (11) 
q = a(u, -0.08)(u^ + 6u* + c) (12) 

The values of a' — 258, a — 182, h — —0.07, c — 0.5, and ut — 0.036 are obtained 
from fitting. The cubic equation is clearly a better fit to the data than the quadratic 
one. 

Fig.® shows the behavior of the grain velocity distribution for the parabolic 
profile. In this case, the saltation regime starts at = 1, below which all particles 
are in creeping motion. Beyond — 1.5, all particles are in the saltation regime. 

4 Discussion 

The dynamics of the granular systems under shear is usually studied for slow or 
quasi static fiow [H] . Slow shearing can lead to the emergence of internal ordering 
of the granular systems. In terms of our study, this means the creeping regime. In 
fact we have also observed size segregation for the case of the logarithmic profile 
when the shear is small enough. Two granular materials differing by their size, 
exhibit a propensity for segregation. Whenever this mixture undergoes a vibration 
or a shearing action, the components tend to separate partially or completely [T^ . 
Here we considered a system with two different sizes, differing in ratios by a factor 2, 
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Figure 6: Mean flux as a function of for a parabolic profile. The solid line is our 
data, and the dashed line is fitting a cubic function to data. Beyond ~ 2.5 the 
effect of the 'top boundary' is revealed, ut — 0.08 is estimated from the simulations. 




Figure 7: Mean flux as a function of u^, near the threshold value ut = 0.08 for a 
parabolic profile. For < 1 we fitted the data with two Eqs.fTTjl , (fT2|) . The solid 
line is the fit using Eq. in]) yielding a = 182, b = -0.07, and c = 0.5. The dashed 
line is the fit using Eq. (fTT|) that gives a' — 258, and ut = 0.036. 
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Figure 8: The behavior of the grain velocity probabihty function for different values 
of for a parabolic profile. 



initially distributed randomly. When this system is affected by the drag force of the 
fluid, the two components of the system separate so that large particles will come 
to the top. In the creeping regime due to the higher collision rate, this separation 
occurs very soon. 

Overall, the present study considers the intermediate region, between the quasi 
static flow and the rapid flow. Due to finite size effects, we use the term 'saltation', 
to mean a motion that takes a particle on a trajectory, however small, raises it from 
the bed of grains. 

We interpret the results using the grain velocity distribution function, P{v). 
With increasing shear we expect to have more and more particles in motion. The 
distribution function, P{v), develops a peak at low velocities, v 2, that gradually 
disappears as u* increases. In case of the logarithmic profile, this peak does not 
broaden. With the diminishing of this low energy peak, a plateau is developed in 
the velocity distribution. This plateau gradually covers a larger range of velocities. 
This plateau corresponds to saltons. The density of particles in the saltation regime, 
for the logarithmic profile, is a monotonically decreasing function of height. We 
attribute this to the comparatively smaller energy input rate; the latter going as 
Qfuj, where qf is the fluid flux. 

For the parabolic profile, the situation is markedly different. We observe that 
as the shear velocity increases, the distribution develops a tail, covering the larger 
velocities. The rate of energy input is much higher in the case of the parabolic 
profile. In the saltation regime (m* > 0.8), it is possible to distinguish a group of 
particles that moves at a higher velocity and height, nearly separating from the rest 
of the particles. In this sense, the density of particles in the saltation regime, is 
qualitatively different from that of the logarithmic velocity profile. This explains 
the sudden velocity spread at > 0.8. 



5 Conclusion 

In this study, we have presented a simple model for sediment transport in creeping 
and saltating motion that produces steady sediment transport. We investigated the 
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steady state fluxes for the parabolic and logarithmic fluid velocity profiles. 

For the logarithmic profile we compared our results with some previous studies. 
Increasing the shear velocity from threshold value Ut, sediment transport sets 
in, first creeping and for larger h*, saltating. In our system, the saltating particles 
can rise up to several times their diameter, and similar to the Ungar and Haff we 
find that the steady state flux increases quadratically with shear velocity. For the 
parabolic distribution the flux rises cubically with increasing the shear velocity. 

The velocity probability distribution function is a good measure to estimate the 
onset of the saltation. It would be interesting to study a system large enough so that 
the grain velocity distribution function can be measured as a function of height. 
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